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Abstract. Morphogens are proteins, often produced in a localised region, whose 
concentrations spatially demarcate regions of differing gene expression in developing 
embryos. The boundaries of expression must be set accurately and in proportion to 
the size of the one-dimensional developing field; this cannot be accomplished by a 
single gradient. Here, we show how a pair of morphogens produced at opposite ends 
of a developing field can solve the pattern-scaling problem. In the most promising 
scenario, the morphogens effectively interact according to the annihilation reaction 
A + B — > and the switch occurs according to the absolute concentration of A or B. 
In this case embryonic markers across the entire developing field scale approximately 
with system size; this cannot be achieved with a pair of non-interacting gradients 
that combinatorially regulate downstream genes. This scaling occurs in a window of 
developing-field sizes centred at a few times the morphogen decay length. 
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1. Introduction 

Morphogen gradients play a crucial role in establishing patterns of gene expression 
during development. These patterns then go on to determine the complex three- 
dimensional morphology that is needed for organism functionality. Because not all 
environmental variation can be controlled, gene patterning must be robust to a variety 
of perturbations, i.e. must compensate for the unpredictable [T]. 

One aspect of this robustness concerns the notion of size scaling Typically, 
gene patterns are established in proportion to the (variable) size of the nascent embryo. 
A dramatic demonstration of this was made recently in the case of Drosophila where 
the posterior boundary of the hunchback gene expression domain was shown to scale 
(to within 5%) with embryo size In the standard model of pattern formation 
in developmental biology cells acquire their positional information by measuring the 
concentration of a morphogen gradient and comparing to some hard-wired set of 
thresholds [3J El El • As the simplest single-source diffusing morphogen gradient with 
fixed thresholds clearly does not exhibit this type of proportionality, it is clear that 
more sophisticated dynamics must be responsible for the observed structures [Zj. 
Unfortunately, little to nothing is known experimentally about how this pattern scaling 
comes about. 

As a first step in deciphering what these more complex processes might entail, we 
study here the issue of how two morphogen gradients, directed from opposite ends of 
a developing field, may solve the pattern-scaling problem Operationally, opposing 
gradients may arise in developing systems in at least two ways. First mRNA, from which 
protein is translated, may be anchored at opposite ends of the region in question. As an 
example, in the Drosophila syncytium an anterior-to-posterior gradient is established 
by the localisation of bicoid mRNA to the anterior, while nanos mRNA localised at the 
posterior defines a reciprocal gradient |B 3 . Second, protein may be secreted by clusters 
of cells, one cluster located at opposite ends of the developing field In both cases we 
will assume that there is some flux of morphogens entering a specific region and assume 
that there is no production in the bulk of the region. We further assume that the 
morphogen reaches neighbouring cells by an effective diffusion process thereby creating 
a gradient Finally, although time-dependent effects in development patterning 

might be important in some contexts [TT], we assume here that a steady-state analysis 
is sufficient insofar as scaling of patterns with system size is concerned. 

We consider two mechanisms in which a pair of morphogen gradients transmits 
size information to the developmental pattern. The first mechanism, which uses the 
concentrations of both gradients combinatorially, is an alternative to the simple gradient 
mechanism • in this mechanism there exist overlapping DNA-binding sites of species 
A and B in the cis-regulatory modules of the target genes. We note that in the 
Drosophila syncyctium some kriippel binding sites overlap extensively with bicoid sites 
|12| IT3] . In this scenario one of the morphogens acts as a transcription factor and the 
other acts as an effective repressor by occluding the binding site of the first. Hence the 
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target gene is switched on according to the relative concentrations of the two species 
[Tlj . In the second mechanism protein B irreversibly inhibits the activity of transcription 
factor A by either directly degrading it or by irreversibly binding to it. The interaction 
is described by the annihilation reaction A + B — > 0. Here we assume that the target 
gene measures the absolute value of the A concentration as in the standard model of 
developmental patterning; the B gradient serves only to provide size information to the 
A concentration field. 

The goal of this work is to study these two possibilities and see the extent to which 
they do in fact solve the pattern scaling problem. To this end we measure the range of 
variables over which scaling is approximately valid. We begin in §Elby pointing out that 
a single gradient in a finite system cannot set markers proportionately in the developing 
field. In § 0] we study the case of two gradients whose binding sites overlap and show 
that approximate scaling then occurs in a restricted fraction of the developing field 
typically located midway between the sources. We then turn to the annihilation model 
of two gradients in § 0] and show that its scaling performance is excellent throughout 
the developing field. 



2. Single Gradient 

Let A be the concentration of the morphogen which in the simplest model obeys 

= D a d 2 x A-(3 a A (1) 

at steady state. Here D a is the diffusion constant of protein A and (3 a is the degradation 
rate. Molecules of A are injected at the left boundary with rate T a and are confined to 
the interval [0, L] by a zero-flux boundary condition 

-D a d x A(0) = T a - -D a d x A(L) = 0. (2) 

The obvious solution is 
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cosh ( — — J = A(L) cosh ( — — ) . (3) 



The length scale A a is defined by A a = JD a /(3 a . 

Let us assume that the boundary between different gene expression regions is 
determined by the position x t at which A equals some threshold value A t . Inverting, 
the expression for the threshold position is 

x t (L) = L - X a cosh" 1 (A t /A(L)) . (4) 

Note that there is a minimum system size for a specific threshold, 

L m = Aasinh" 1 j , (5) 

such that Xt(L m ) = L m . When L — x t ^> A a the concentration profile becomes purely 
exponential and where 
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Clearly, the function x t (L) starts out at L m (which is greater than Xoo), is monotonically 
decreasing, but is bounded below by x^; in other words x t is always greater than 
as the effect of the zero-flux boundary condition is to make xt larger than it would be 
in the absence of the boundary. Fig. U shows the variation of x t /L with L for three 
different values of the threshold concentration A t . There is no extremum where this 
ratio becomes locally L- independent. 



3. Combinatorial model 



We next ask whether a molecular mechanism, which compares the concentrations of 
two gradients rather than reading the absolute value of one or both of them, can lead 
to gene expression boundaries which scale with system size. We consider two opposing 
gradients A and B described by 

= D a d 2 x A - (3 a A (7) 
= D b d 2 x B - (3 b B (8) 

at steady state. The boundary conditions are 

-D a d x A(0) = r o ; -D a d x A{L) = 

-D b d x B(0) = 0; —D b d x B(L) = -T b . [ ' 

In this scenario, the gene expression boundary will be determined by a critical 
concentration ratio r which occurs at the position x r defined by A(x r ) = rB(x r ). 

Just as in the one-gradient case, one can distinguish between relatively small 
systems (for which the no-flux boundary conditions matter) and large systems, 
depending on how big L is compared to the decay lengths A«. For sufficiently large 
L, the gradients of A and B are purely exponential, A = A(0) exp(— x/X a ) and 
B = B(L) exp(— (L — x)/X b ), and the gene expression boundary is given by 

Xr = ( L ± L c) ( 10 ) 



where 



A,, 



In 



'rB(L) 



A{0) . 

Consider a gene whose cis-regulatory module contains overlapping A and B binding 
sites. This gene will have a particular threshold ratio r and a concomitant value L c (r). 
Then for sufficiently large developing fields L ^> L c (r) the combinatorial mechanism 
sets the boundary of expression of the gene at the relative location 

T - < 12 > 

in a size- invariant manner. This position is also insensitive to source-level fluctuations, 
which only enter in L c .. In a system in which the degradation lengths of the two 
morphogen gradients are comparable, x r /L will be close to 1/2. 
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Although this model can achieve some degree of size-scaling near the centre of the 
developing field, from Eq. (fTU|) it is clear that the variation of x r /L with L increases 
as x r /L deviates from the aforementioned asymptotic value. This will happen either as 
the size of the system is made smaller, or even at fixed L if we try to make the threshold 
point x r approach the edges of the developing field. We have in mind a situation where 
multiple genes need to be regulated, each at different points along the developing field; 
each gene will have its own value of r and hence its own value of x r . In the previous 
limit, there is no variation in x r with r and this cannot be accomplished; therefore we 
need to rely on finite L c / L effects. To proceed, we must more carefully characterize the 
variation of x r /L with L for all positions in the developing field. Since Eq. (fTT)|) becomes 
inaccurate close to the edges of the developing field we return to the expression for A 
in Eq. (and a similar one for B) and obtain the following implicit equation for x r 

ML) = rB(0) 
cosh(x r /A a ) cosh((L — x r )/\ b ) 

valid for a finite system. It will be critical to identify what happens to x r when the 
length L is made smaller. Notice that there is a different behavior depending on which 
of A(L) and rB(0) is larger. Specifically, if A(L) is larger, there will be a smallest length 
below which x r given by this formula becomes larger than L; this length is given by 



L*(r) = \ a cosh" 



-i ( ML) \ 
\rB(0) 



If, on the other hand, the ordering is reversed, below the length scale 

we obtain negative values for x r . Representative x r /L curves are shown in Fig. |2]for 
the case of equal decay lengths X a = X b . 

Consider now a developing field of size L subject to a natural variation in size of 
L ± pL with < p < 1 . The variation in the fractional position at which a gene is 
turned on is then given by 

fXr\ = x r (L-pL) _ x r (L+pL) 

\LJ ~ L-pL ~ L + pL 1 ' 

We show in Fig. Efa), again for the equal decay length case, the dependence of S (x r /L) 
on normalised position x r /L in the developing field for L = 4. As expected the variation 
is largest (in magnitude) at the boundaries and vanishes at that position x ro for which 
the critical length L c {tq) vanishes. Defining an arbitrary scaling criterion according to 

5{x/L)<5%, (15) 

one sees that the combinatorial model achieves scaling only in the central region of the 
developing field between about 30% and 70% of L. 

Near the edges of the developing field the variation 5(x r /L) is about 14%. Since 
the slopes of the x r /L curves at x r /L = 1 become flatter as L is increased (see Fig. EJ), 
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one might wonder whether operating at larger system sizes will decrease this variation. 
However at larger system sizes the flattening effect is offset by the fact that one must 
sample larger and larger portions of the x r /L curve. The extent to which these effects 
cancel is shown in Fig. 0(b) where we show the variation 5 (x r /L) closest to the right 
boundary of the developing field as a function of L. The variation decreases with L, 
but an elementary calculation reveals that it has the lower bound p/{l + p). For a 
percentage variation p = 10% in system size this lower bound is about 9%. We conclude 
that increasing system size is not sufficient to make the combinatorial model, with 
A a = Xb, meet the scaling criterion throughout the developing field. 

A further difficulty with the combinatorial model is its susceptibility to small- 
molecule-number fluctuations. In general, we must expect L c of order A, since we 
cannot independently adjust the morphogen sources for the multiple genes that need to 
be controlled. In fact, the natural interpretation of r as being due to binding differences 
between different transcription factors suggest that L c would vary significantly. In such 
cases the limit L ^> L c (tq) would force the comparison point x r far down the profile from 
the source; having enough molecules at this point to affect the necessary DNA binding 
would then place a severe constraint on source strengths. In this regard a combinatorial 
mode of action may favour power-law (resulting e.g. from nonlinear degradation [To] ) 
over exponential profiles as the former have greater range than the latter, but this 
remains to be studied. 

4. Annihilation model 

We return to the standard model of morphogenesis in which cell-fate boundaries are 
determined according to the position at which a single morphogen crosses a threshold 
concentration. We couple this gradient to an auxilary gradient directed from the 
opposite end of the developing field. We then ask under what conditions the primary 
gradient may scale with system size. 

We consider two species of morphogen, A and B, in a one- dimensional system 
of length L with As and Bs injected at opposite ends of the system. The boundary 
conditions are as in § El The species interact according to the annihilation reaction 
A + B —>■ 0. In a mean-field description the kinetics is described by the reaction-diffusion 
equations 



where k is the annihilation rate constant. Later, we will consider more complex models 
which incorporate non-linear degradation or non-linear (i.e. concentration-dependent) 
diffusion. 

This system of equations, with fluxes r a = Ft, = T and without any decay, 
was considered by Ben-Nairn and Redner |T^j. They determined the steady-state 
spatial distribution of the reactants and of the annihilation zone which they chose 



d t A = D a d 2 x A - (3 a A - kAB 
d t B = D b d 2 x B - (3 b B - kAB 



(16) 
(17) 
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to be centred in the interval [0, L\. The annihilation zone is roughly the support of 
R(x) = kA(x)B(x) or, put another way, that region where the concentration of both 
species is appreciable. With the aid of a rate-balance argument, they showed that the 



is proportional to T 2 / 3 when w <^ L. 

Our goal is to understand the relation of the steady-state concentration profiles to 
the system length L. It is convenient to identify the point x e in the annihilation zone 
where the profiles cross, A(x e ) = B(x e ). In the original Ben-Nairn — Redner model, the 
reaction-diffusion equations yield no unique value for x e ; instead x e can lie anywhere in 
the interval [0, L] depending on the choice of initial condition. To see this consider the 
following rate-balance argument. Since the particles annihilate in a one-to-one fashion 
the flux of each species into the annihilation zone must be equal. But this condition 
does not determine x e uniquely because these fluxes are always equal to the input fluxes 
at the boundaries. Similarly, the model without degradation cannot support steady 
states with unequal boundary fluxes. If, however, we now add degradation terms to the 
steady-state equations, then the flux of each species into the annihilation zone is the 
flux into the system less the number of degradation events that happen before reaching 
the zone. Thus, the flux of each species into the annihilation zone now depends on the 
location x e and so there is only one value of x e which balances the fluxes. As we will 
see, our models will always contain unique steady-state solutions. 

A rough estimate of the concentration in the annihilation zone and of the width of 
the zone can be obtained using the original Ben-Naim-Redner rate-balance argument 
[To] . We identify three spatial regions: the first where A is in the majority; the second the 
annihilation zone; and the third where A is in the minority. Assume the concentration 
of As in this latter region is negligible compared with that in the other two regions. The 
concentration of As in the annihilation zone should then be of the order of the slope of 
the concentration profile in the annihilation zone times the width w. The slope of the A 
profile in this region is proportional to j e /D a , where j e is the equal flux of As or Bs into 
the annihilation zone. Therefore the concentration in the annihilation zone A e = A(x e ) 
is 



If we ignore the loss of A particles in the annihilation zone (valid for small w) , then the 
number of annihilation events per unit time kA 2 e w should equal the flux j e . Balancing 
these two rates gives j e ~ k(j e w/D a ) 2 w. Hence the width of the annihilation zone scales 
as 



In what follows, we will be mostly interested in taking k large enough to give a very 
small w. 



width w of the annihilation zone scales as T and that the concentration in this zone 



A e ~ j e w/D a . 





(19) 
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5. The high-annihilation-rate limit 

We now explicitly assume that the parameters lie in the limit where w -C min{a; e , L—x e }. 
This limit has the considerable advantage that the A-B system may be decoupled by 
replacing the coupling term kAB by a zero-concentration boundary condition at x e . In 
this approximation the concentration of the A subsystem satisfies 

= D a d 2 x A - (3 a A (20) 

subject to the boundary conditions — D a d x A(0) = T a and A(x e ) = . The solution to 
this equation is 

A(x) = - X f. a m \ smh (^f^) = A sinh (^f^) (21) 



D a cosh(x e /A a ) \ A a J V A a 



where as before A a = y D a / f3 a . A* is a characteristic concentration of the A field related 
to the slope of the A field at x e according to A* = —X a d x A(x e ). The flux of A particles 
is 

ja(x) = j a {x e ) COsh ~) ' ^ 

where the flux into the annihilation zone j a (x e ) is given by j a {x e ) = T a / cosh(x e /A a ). 
Substituting this into Eq. yields the scaling function of the annihilation zone width 
for the case of linear degradation 

w ~ w [cosh(a; e / A a )] 1/3 . (23) 

1 /3 

Here Wq ~ (D 2 /T a k) is the width of the annihilation zone in the absence of 
degradation [TH]. Note that we may also substitute this expression for j a (x e ) into 
Eq. (|18jl obtaining A e ~ w/ cosh(x e /A a ). One can then verify that A e is much smaller 
than A(0) whenever w<Ci e and hence approximating this as a zero boundary condition 
is self-consistently valid. 

The .B-subsystem can be treated similarly, except that the length of the subsystem 
in this case is L — x e . The only dependence on the annihilation rate k in the inequality 
t«Ci e occurs in wq. Hence this limit is equivalent to the high-annihilation-rate limit 
k ^> k , where the threshold value k of the annihilation rate is given by 
D\ cosh(q; e /A a ) 

h ~ WW ■ ( ' 

We determine the annihilation zone location by balancing fluxes into the zone, 
3a{ x e) = —jb(L — x e ). This leads to the following equation for x e 
V V 

(25) 



cosh(f^) cosh(^f)' 

In the special case A a = A;, this equation coincides with the implicit definition of x r 
(with r = 1) which arose in the combinatorial model (see Eq. (|13|1). As in that model 
there is a smallest length L* defined by 

T„ 



Embryonic Pattern Scaling 



9 



if T a > T b and by 

L^A^cosh- 1 

if the flux ordering is reversed. As our entire treatment of the annihilation zone 
only makes sense if < x e < L, we must always choose L > L*. A comparison of 
the numerical solution of the full model with the results of the large-annihilation-rate 
approximation is shown in Fig. 

Once we know x e (L) and A(x), we can proceed to determine the qualitative features 
of the x t (L) function with a view to identifying the region of system sizes where x t ~ L. 
Inverting Eq. (|21|) we find 

Xt = x e — X a sinh -1 77 (26) 

where 

V = A/A,. (27) 

Note that x t depends on L only through its dependence on x e and the function x t (x e ) is 
monotonically increasing. Obviously x t < x e . In the limit of sufficiently large x e , we can 
replace the inverse hyperbolic function with a logarithm and obtain the simpler form 

x t &x e - A a lm». (28) 

Here, n « ^e Xe ^ Xa , and x t approaches its asymptotic value Xoo ~ A a ln(yl(0) / 'A t ) from 
below. This is of course the answer one would obtain in the absence of any auxiliary 
gradient. 

Now, imagine reducing L and hence x e from its just-mentioned asymptotic regime 
and plotting the ratio x t /L. For the case T a > Tb, x e will eventually hit L followed 
shortly thereafter by Xt/L hitting unity. There is no reason why this curve should 
exhibit a maximum, and a direct numerical calculation for k = 100 (shown in Fig. [5J) 
verifies this assertion. The situation is dramatically different, however, for the case 
of Tfe > r o . Now x e must approach zero, implying that at some larger L we have 
Xt=0. The curve x%jL now exhibits a maximum, as is again verified by direct numerical 
calculations using both the large-annihilation-rate approximation and also just solving 
the initial model with no approximations whatsoever (see Fig. EJ). Near the peak of the 
curve we have scaling with system size. For completeness, we also present in Fig. [7|the 
results for equal fluxes. 

To compare the scaling performance of the annihilation model with that of the 
combinatorial model we show in Fig. [Sj the dependence of the variation 5(x t /L) on 
normalised position Xt/L in the developing field for L = 4. One sees that, according to 
our scaling criterion in Eq. (JT5)l , the annihilation mechansim can easily set markers scale- 
invariantly throughout a developing field whose size is a few decay lengths. Furthermore 
at such system sizes a range of threshold values spanning two orders of magnitude 
(A t = 0.01 — 0.7) is sufficient to cover the entire developing field (see Fig. EJ). Such 
a modest variation in concentration makes the annihilation model less susceptible to 
small-molecule-number fluctuations than the combinatorial model. 
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6. Discussion 

We have considered two scenarios in which a pair of oppositely directed morphogen 
gradients are used to set embryonic markers in a size-invariant manner. In the simplest 
scenario, in which the gradients interact only indirectly through overlapping DNA- 
binding sites, exponentially distributed fields achieve perfect size scaling at a normalised 
position A a /(A a + A&) determined only by the morphogen decay lengths A a and A&. For 
equal decay lengths, the accuracy with which this model can set markers size-invariantly 
decreases as the boundaries of the developing field are approached. At the boundaries 
the accuracy can be no better than 5(x r /L) = p/(l + p) where p is the percentage 
variation of the field size. In the second model A and B are coupled via the reaction 
A + B — > and the embryonic markers are set by a single gradient with the second 
gradient serving only to provide size information to the first. In this scenario, it is easy 
to arrange parameters such that scaling occurs with an accuracy 8{x t /L) better than 
5% over the entire developing field for field sizes of only a few decay lengths. 

In practice a given morphogen may play both roles in patterning, setting markers 
in a strictly concentration-dependent manner at some locations in the developing field 
and in a combinatorial fashion at other locations [I2|. The annihilation model naturally 
sets markers via the gradient whose source is closest to the marker ^7j, whereas the 
combinatorial model is better suited to setting markers in the vicinity of the midpoint 
of the developing field where the variation 5(x r /L) is smallest. As the variation 5(x/L) 
has a qualitatively different dependence on x/L in either measurement of this 

curve in a developmental system may distinguish between the mechanisms. 

The origin of the scaling form f(x/L) which arises in the strong-coupling limit 
of the annihilation model is the effective boundary condition A(x e ) = 0. In the case 
T b > T a (see Fig- IHI) the x t /L curve has a maximum because at small L (L ~ L*) it tends 
to zero along with x e /L while at large L [L L*) it is bounded above by Xoq/L. In the 
k <C k limit, on the other hand, the zero-concentration effective boundary condition is 
replaced by a zero- flux boundary condition j a {L) = which can never induce the Xt ~ L 
scaling. 

This approach makes it clear why the scaling occurs at intermediate values of 
L. Once we reach the non-overlapping limit where the two fields do not effectively 
communicate, the threshold is set by the A profile alone; we have already seen that 
this cannot give any scaling. For L too small, the annihilation-zone width w becomes 
comparable to x e , there is no effective boundary condition and again scaling fails. In 
fact, if one looks at the expression for w/x e , namely 

W_ Wq ( COsh(x e /A a ) \ 1/3 

Xe ~ A a ^ (xJKf J ' 1 j 

(where wq ~ {D^/Y a k) 1 ^) one sees that the maximum in Xt/L occurs close to the 
minimum of w/x e which is reached at x e /X a ~ 3. 

So far we have used linear degradation and simple diffusion in the annihilation 
model. However, it should be clear from the above arguments that in fact this mechanism 
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is rather robust to changing the nature of the individual gradients. For example, let 
us consider quadratic degradation. In the limit that the system size is so big as to 
render the coupling term kAB irrelevant the A and B profiles reduce to power laws, 
A = a/(x + e a ) 2 and B = b/(L — x + e^) 2 . The corresponding L-independent threshold 
position Xoo is given by 



An argument, similar to one presented earlier for linear degradation, reveals the fact 
that x e will be forced to zero for sufficiently small L if 1^ > T a ; this indicates again that 
to the extent we can believe the large-annihilation-rate approximation, there will be a 
maximum in the x t /L curve. This is illustrated for one specific choice of parameters 
in Fig. EKa) • The maximum again takes place roughly where L becomes so small as 
to cause the annihilation- zone width to approach x e . Repeating the derivation of w 
outlined in § El but using a power law instead of hyperbolic sine we obtain 



This expression is a good qualitative description of the exact w/x e shown in Fig. Efa) 
and diverges when L — > as in the case of linear degradation. Notice that scaling is lost 
when w — > x e even though the rate of the annihilation reaction becomes large (Fig. 0(b)). 
Finally, one can also ask about the effect of making the diffusion constant concentration 
dependent. This type of effect can arise whenever the morphogen reversibly binds to 
buffers that differ in mobility from the pure molecule. Fig. El illustrates the behavior 
under the simplest assumption, namely that the diffusion constant varies linearly with 
concentration for both the A and B fields. Aside from sharpening the transition from 
the asymptotic non-interacting regime to the regime where x e approaches zero (as L is 
lowered), the basic phenomenology is unchanged. 

The focus of our work has been the scaling issue. However, we should not lose 
track of the other requirement for developmental dynamics, namely that the system 
be relatively robust to fluctuations in parameters such as source fluxes. Fig. ITTTa) 
presents data regarding the variation of x t with r a and in the annihilation model. 
For simplicity the data is presented for the case of equal decay lengths, A a = A& = A. 
The basic conclusion is that the coefficient of variation Xi, defined as 



starts at 1/2 at A t — and then asymptotes to either 1 for variations in T a or 
zero for variations in Tb- These asymptotic values are of course precisely the results 
obtained for the one-exponential-gradient model. The fact that the Xi a ^ small x t is 
1/2 can be understood by noting that in this limit xt is just x e , which can easily be 
shown to be approximately (i.e. for large enough L) given by x e ~ 0.5(L ± L c ) with 




(30) 




(31) 




(32) 
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L c = A |ln(r&/r o )|. With this approximation for x e and taking differentials of x t we 



obtain 



Xa 



Xb 



V 



1 

2 + v^T^ 



1 - 



7/ 



2 

tanh 



- 1 anli ( 



(33) 

t m , ; (34) 

V 1 + rf 

where, as before, r] = A t /A*. These are good approximations at all values of 77 for 
percentage variations in source fluxes as large as 5% (see Fig. ITTT a)). The reduction of 
the x values from unity represent an increase in system robustness as compared with 
the single-exponential-gradient model, albeit with a new sensitivity to the B gradient. 
For comparison we also show in Fig. fTTT b) the coefficient of variation which arises in 
the single-gradient model. The approximation to Xa in this case is given by 

V 



Xa 



1 



(35) 



where now 77 is defined by 77 = A t /A(L). Notice that the effect of the boundary (77 J, 1) 
is to increase the sensitivity of the gradient to variations in the source flux over that for 
a simple exponential. 



7. Conclusions 



In this paper we have shown that coupling two oppositely directed morphogen gradients 
allows patterns to be set in approximate proportion to the size of the developing field. 
We have considered two coupling mechanisms, the most effective of which couples 
the gradients via a phenomenological annihilation reaction. Such a mechanism can 
set boundaries of gene expression across the developing field with a small sample- 
to-sample variation in the normalised position of the boundaries. In this scenario, 
there is no magic bullet which ensures either exact scaling or complete robustness. 
Instead, the effective boundary condition created by the annihilation reaction allows 
for approximate scale invariance to emerge in one reasonably-sized range of parameter 
space and similarly lowers the sensitivity of any threshold to source-level fluctuations. 
Presumably, one could obtain even more robustness and scaling, and possibly even 
temperature compensation (see for example Ref. P-8J), via the introduction of yet 
additional interactions. 

After completion of this work we became aware of similar work in which the 
annihilation model was applied to pattern scaling in the early Drosophila embryo [TTH l2U|. 
In contrast to the numerical analysis carried out by the authors of Ref. |2Hj for the specific 
case of the bicoid morphogen, we have presented here a more general analytic framework 
which allows for a natural explanation of pattern scaling at intermediate developing-field 
sizes and of filtration of source-level fluctuations. In particular our work provides an 
explanation for the fact that they found pattern scaling at L approximately 4-5 times 
the decay length A. We note that recent work showing that only bicoid binding sites 
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are needed for scaling provides further support for an annihilation mechanism in the 
bicoid-hunchback problem [21] . 
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Figure 1. Dependence of normalised xt on system length L for the case of a single 
gradient. The solid lines are the analytic expressions Eq. J3J for x t /L for values of 
the threshold concentration equal to (from top to bottom) At — 0.01,0.1,0.7. Dashed 
lines are Xoo/L curves as given by Eq. @. All parameters are unity unless otherwise 
stated. 
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Figure 2. Dependence of normalised x r on system length L in the combinatorial two- 
gradient model, as given by Eq. for values of the threshold ratio equal to (from 
top to bottom) r = 1CP 5 , 1CP 3 , 10 _1 ,0.5, 1. All parameters are unity unless otherwise 
stated. 
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Figure 3. (a) Dependence of the variation 8{x r /L) on normalised position x r /L in 
the developing field for L — 4 and a percentage change in system size of 10%. (b) The 
variation S(x r /L) closest to the right boundary of the developing field as a function 
of L (solid line). At each L we have chosen the target gene whose threshold ratio r 
satisfies L*(r) = L — pL. The variation in the fractional position at which this gene 
is turned on is then given by6(^) = l— ^jjq^f^ ■ The dashed line is an asymptotic 
expression. The horizontal (red) line is the limiting value p/(l+ p) of the solid and 
dashed curves. 
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Figure 4. (a) Comparison of the high-annihilation-rate approximation (Eq. I)21[l. solid 
line) with the numerical solution (circles) of the full annihilation model (Eqs. 
and l|17ll ). The annihilation zone is the reaction front R(x) — kA(x)B(x). All 
parameters are unity unless otherwise stated, (b) A(x) plotted on a logarithmic scale 
in the cases k = 0.01 and k = 100. Note the crossover from slow decay in the A-rich 
region to fast decay in the B-rich region in the case k = 100. All parameters are unity 
unless otherwise stated. 




Figure 5. Dependence of normalised Xt on system length L with k = 100 and Ft, < T a . 
The plus signs, circles and crosses are numerical solutions of Eqs. Ijl6(l and (|1 T|) for 
values of the threshold concentration equal to (from top to bottom) A t = 0.01, 0.1, 0.7. 
The solid lines are the corresponding analytic expressions Eq. (|26J) obtained in the 
high- annihilation-rate limit. Dashed lines are x^jL curves as given by Eq. ©. All 
parameters are unity unless otherwise stated. 
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Figure 6. Dependence of normalised Xt on system length L with k = 100 and Ft, > T a . 
The plus signs, circles and crosses are numerical solutions of Eqs. Ijl6(l and (|1 T|) for 
values of the threshold concentration equal to (from top to bottom) A t = 0.01, 0.1, 0.7. 
The solid lines are the corresponding analytic expressions Eq. (|26J) obtained in the 
high- annihilation-rate limit. Dashed lines are x^jL curves as given by Eq. ©. All 
parameters are unity unless otherwise stated. 
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Figure 7. Dependence of normalised Xt on system length L with k = 100 and Ft, = T a . 
The plus signs, circles and crosses are numerical solutions of Eqs. Ijl6(l and (|1 T|) for 
values of the threshold concentration equal to (from top to bottom) A t = 0.01, 0.1, 0.7. 
The solid lines are the corresponding analytic expressions Eq. (|26J) obtained in the 
high- annihilation-rate limit. Dashed lines are x^jL curves as given by Eq. ©. All 
parameters are unity unless otherwise stated. 
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Figure 8. The dependence of the variation 5(xt/L) on normalised position Xt/L in 
the high-annihilation-rate approximation of the annihilation model. Positions to the 
left of x e are set by the A gradient while positions to the right are set by the B gradient. 
All parameters are unity unless otherwise stated. 
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Figure 9. (a) Dependence of normalised xt on system length L in the case of quadratic 
degradation. The plus signs, circles and crosses are numerical solutions of the full 
annihilation model for values of the threshold concentration equal to (from top to 
bottom) A t = 0.01,0.1,0.7. The dashed lines are x x /L curves as given by Eq. J3UJI. 
Also shown (cyan diamonds) is the ratio of the full width at half maximum w to 
the comparison point x e . (b) The dependence of the amplitude R ma x of the local 
annihilation rate R(x) = kA(x)B(x) on system length L. In (a) and (b) all parameters 
are unity unless otherwise stated. 
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Figure 10. Dependence of scaled threshold position x t /L on system length L in 
the simplest case of nonlinear diffusion, D a = S a A and £>{, = SbB. The degradation 
terms are linear. The plus signs, circles and crosses are numerical solutions of the 
full annihilation model for values of the threshold concentration equal to (from top to 
bottom) At = 0.01, 0.1, 0.7. Dashed lines are corresponding curves in the case k = 0.01. 
All parameters are unity unless otherwise stated. 
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Figure 11. Sensitivity of the threshold position xt to infinitesimal variations in the 
source fluxes T a and Tj, in (a) the annihilation model (constant diffusion constant and 
linear degradation) and (b) the single-gradient model. The coefficient of variation Xi 
is defined by Eq. llol'l) in the text. The data plotted as plus signs and circles were 
obtained by solving numerically the full model, while the dashed lines represent (from 
top to bottom) Eqs. I|34|) and i|35|) . All parameters are unity unless otherwise 

stated. 



